library(librarian)
shelf(tximport, tidybulk, tidySummarizedExperiment, ggpubr, tidyverse)

wu21meta <- read_tsv('Archive/covid19/data/wu-chen-samples') |>
  select(id1, clinic)

t <- read_rds('~/append-ssd/nextflowing/rnaseq_wu2021_batch1/salmon/salmon.merged.gene_counts_length_scaled.rds')

t2 <- read_rds('~/append-ssd/nextflowing/rnaseq_wu2021_batch2/salmon/salmon.merged.gene_counts_length_scaled.rds')

t3 <- read_rds('~/append-ssd/nextflowing/rnaseq_wu2021_batch3/salmon/salmon.merged.gene_counts_length_scaled.rds')

tdb <- tidybulk(t)

bk_list <- list(t,t2,t3) |>
  map(\(x)x |>
        tidybulk() |>
  mutate(.feature = gene_name) |>
  aggregate_duplicates(.transcript = gene_name) |>
  select(1:4)
)

tdb <- bk_list |>
  list_rbind()

tdb <- tdb |>
  mutate(id1 = str_extract(.sample, 'P...')) |>
  left_join(wu21meta)

tdb |>
  filter(is.na(clinic)) |>
  dplyr::count(id1)

tdb <- tdb |>
  filter(!is.na(clinic))

seg_anno <-
'x,x2,y
1,2,12
1,3,13
1,4,14' |> read_delim()

pval_anno <-
  'x,y,p
1.5,12.3,*
2,13.3,*
2.5,14.3,***' |> read_delim()

tdb |>
  filter(.feature == 'HMCES') |>
  ggplot(aes(clinic, abundance, color = clinic)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  scale_color_brewer(type = 'div', direction = -1) +
  geom_segment(data = seg_anno, aes(x = x, xend = x2, y = y, yend = y), color = 'black') +
  geom_text(data = pval_anno, aes(x = x, y = y, label = p), color = 'black', size = 5) +
  theme_pubr() +
  labs(x = 'Severity', color = 'Severity', y = 'Normalized expression',
       title = 'HMCES expression in PBMC',
       subtitle = 'Wu 2021 ERP127339')

tdb <- tdb |>
  tidybulk(.sample = .sample,
           .transcript = .feature,
           .abundance = counts,
           .abundance_scaled = abundance) |>
  identify_abundant(factor_of_interest = clinic)

tdb |>
  ggplot(aes(abundance, color=clinic)) +
  geom_density() +
  scale_x_log10()

pca_res <- tdb |>
  filter(!is.na(abundance)) |>
  reduce_dimensions(method="PCA", .dims = 2) |>
  pivot_sample()

pca_res |>
  ggplot(aes(PC1, PC2, color = clinic)) +
  geom_point() +
  scale_color_brewer(type = 'div')

tdb |> write_rds('mission/HMCES_HIV/wu21.tidyblk.rds')

tdb <- read_rds('mission/HMCES_HIV/wu21.tidyblk.rds')

tdb$clinic |> unique()

dif_tdb <- tdb |>
  test_differential_abundance(~ 0 + clinic,
                              contrasts = 'clinic1_asymptomatic-clinic3_severe')

edg_res <- dif_tdb |>
  keep_abundant() %>%
  pivot_transcript()

edg_res |>
  filter(.feature == 'HMCES') |> DT::datatable()

